{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MeanHamilMinimizer native scipy \n",
    "\n",
    "The purpose of this notebook is to describe a naive but pedagogical,\n",
    "first baby step, in the implementation\n",
    "of what is called in Qubiter the Mean Hamiltonian Minimization problem.\n",
    "\n",
    "The qc history of this problem started with quantum chemists planning to\n",
    "use on a qc the phase estimation algo invented by Kitaev? (an algo that\n",
    "is also implemented in Qubiter) to estimate the energy levels (\n",
    "eigenvalues) of simple molecules, initially H2. Then a bunch of people\n",
    "realized, heck, rather than trying to estimate the eigenvalues of a\n",
    "Hamiltonian by estimating the phase changes it causes, we can estimate\n",
    "those eigenvalues more efficiently by estimating the mean value of that\n",
    "Hamiltonian as measured empirically on a qc. Basically, just the\n",
    "Rayleigh-Ritz method, one of the oldest tricks in the book. One of the\n",
    "first papers to propose this mean idea is\n",
    "https://arxiv.org/abs/1304.3061 Their algo is commonly referred to by\n",
    "the ungainly name VQE (Variational Quantum Eigensolver) VQE was\n",
    "originally applied to do quantum chemistry with a qc. But now Rigetti\n",
    "and others have renamed it hybrid quantum-classical quantum computing\n",
    "and pointed out that it's an algo that has wide applicability, not just\n",
    "to quantum chemistry.\n",
    "\n",
    "The idea behind hybrid quantum-classical is very simple. One has a\n",
    "classical box CBox and a quantum box QBox. The gates of QBox depend on N\n",
    "gate parameters. QBox sends info to CBox. CBox sends back to QBox N new\n",
    "gate parameters that will lower some cost function. This feedback\n",
    "process between CBox and QBox continues until the cost is minimized. The\n",
    "cost function is the mean value of a Hamiltonian which is estimated\n",
    "empirically from data obtained from the qc which resides inside the QBox.\n",
    "\n",
    "To minimize a function of N continuous parameters, one can use some\n",
    "methods like simulated annealing and Powell that do not require\n",
    "calculating derivatives, or one can use methods that do use derivatives.\n",
    "Another possible separation is between methods that don't care which\n",
    "local minimum they find, as long as they find one of them, and those\n",
    "methods that try to find the best local minimum of them all, the so\n",
    "called global minimum. Yet another separation is between methods that\n",
    "allow constraints and those that don't.\n",
    "\n",
    "Among the methods that do use derivatives, the so called gradient based\n",
    "methods only use the 1st derivative, whereas other methods use both\n",
    "first (Jacobian) and second (Hessian) derivatives. The performance of\n",
    "those that use both 1st and 2nd derivatives degrades quickly as N grows.\n",
    "Besides, calculating 2nd derivatives is very expensive. Hence, methods\n",
    "that use the 2nd derivatives are practically useless in the neural\n",
    "network field where N is usually very large. In that field, gradient\n",
    "based methods rule.\n",
    "\n",
    "A method that uses no derivatives is Powell. A gradient based method\n",
    "that is designed to have a fast convergence rate is the Conjugate\n",
    "Gradient (CG) method. Another gradient based method is back-propagation\n",
    "(BP). BP can be implemented as distributed computing much more easily\n",
    "than other gradient based methods so it is favored by the most popular\n",
    "computer programs for doing distributed AI, such as PyTorch and\n",
    "Tensorflow.\n",
    "\n",
    "Qubiter can perform minimization using various minlibs (minimization \n",
    "software libraries) such as 'scipy', 'autograd', 'pytorch', 'tflow'. It \n",
    "can also use various devices (aka simulators or backends), either \n",
    "virtual or real, to do the minimization. \n",
    "\n",
    "Non-scipy minlibs implement backprop.\n",
    "\n",
    "The 'scipy' minlib is a wrapper for the scipy function\n",
    "`scipy.optimize.minimize`. This scipy umbrella method implements many\n",
    "minimization methods, including Powell and CG.\n",
    "\n",
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html\n",
    "\n",
    "By a native device, we mean one that uses Qubiter native simulators like\n",
    "SEO_simulator.\n",
    "\n",
    "\n",
    "So, without further ado, here is an example of the use of \n",
    "class `MeanHamilMinimizer` with a scipy minlib and native device."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{\\bra}[1]{\\left\\langle{#1}\\right|}$\n",
    "$\\newcommand{\\ket}[1]{\\left|{#1}\\right\\rangle}$\n",
    "test: $\\bra{\\psi}M\\ket{\\phi}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First change the directory to the Qubiter directory and add it to the path environment variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/rrtucci/PycharmProjects/qubiter/qubiter/jupyter_notebooks\n",
      "/home/rrtucci/PycharmProjects/qubiter\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "print(os.getcwd())\n",
    "os.chdir('../../')\n",
    "print(os.getcwd())\n",
    "sys.path.insert(0,os.getcwd())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we construct a simple 4 qubit circuit that depends on two placeholder variables\n",
    "`#1` and `#2`. These are the continuous variables that we will vary\n",
    " to minimize a cost function. The cost function will be specified later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loaded OneQubitGate, WITHOUT autograd.numpy\n"
     ]
    }
   ],
   "source": [
    "from qubiter.adv_applications.MeanHamil_native import *\n",
    "from qubiter.adv_applications.MeanHamilMinimizer import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_qbits = 4\n",
    "file_prefix = 'mean_hamil_native_test'\n",
    "emb = CktEmbedder(num_qbits, num_qbits)\n",
    "wr = SEO_writer(file_prefix, emb)\n",
    "wr.write_Rx(2, rads=np.pi/7)\n",
    "wr.write_Rx(1, rads='#2*.5')\n",
    "wr.write_Rn(3, rads_list=['#1', '-#1*3', '#2'])\n",
    "wr.write_Rx(1, rads='-my_fun#2#1')\n",
    "wr.write_cnot(2, 3)\n",
    "wr.close_files()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code above wrote inside Qubiter's `io_folder`, two files,\n",
    "an English file and a Picture file. We next ask the writer object \n",
    "wr to print those two files for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style='font-family:monospace'><tr><td style='border-right:1px solid red;'>1</td><td style='text-align:left;'><pre>ROTX\t25.714286\tAT\t2</pre></td></tr><td style='border-right:1px solid red;'>2</td><td style='text-align:left;'><pre>ROTX\t#2*.5\tAT\t1</pre></td></tr><td style='border-right:1px solid red;'>3</td><td style='text-align:left;'><pre>ROTN\t#1\t-#1*3\t#2\tAT\t3</pre></td></tr><td style='border-right:1px solid red;'>4</td><td style='text-align:left;'><pre>ROTX\t-my_fun#2#1\tAT\t1</pre></td></tr><td style='border-right:1px solid red;'>5</td><td style='text-align:left;'><pre>SIGX\tAT\t3\tIF\t2T</pre></td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wr.print_eng_file(jup=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style='font-family:monospace'><tr><td style='border-right:1px solid red;'>1</td><td style='text-align:left;'><pre>|   Rx  |   |</pre></td></tr><td style='border-right:1px solid red;'>2</td><td style='text-align:left;'><pre>|   |   Rx  |</pre></td></tr><td style='border-right:1px solid red;'>3</td><td style='text-align:left;'><pre>R   |   |   |</pre></td></tr><td style='border-right:1px solid red;'>4</td><td style='text-align:left;'><pre>|   |   Rx  |</pre></td></tr><td style='border-right:1px solid red;'>5</td><td style='text-align:left;'><pre>X---@   |   |</pre></td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wr.print_pic_file(jup=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The circuit depends on a placeholder function called `my_fun`.\n",
    "This function will remain fixed during the\n",
    "minimization process but it must be defined \n",
    "and it must be passed in inside an input dictionary to \n",
    "the constructor of class `MeanHamil`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_fun(x, y):\n",
    "    return x + .5*y\n",
    "\n",
    "fun_name_to_fun = {'my_fun': my_fun}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One must also pass into the constructor of `MeanHamil`, the value of a Hamiltonian called `hamil`.\n",
    "Qubiter stores `hamil` in an object of the class `QubitOperator`\n",
    "from the open source Python library `OpenFermion`.\n",
    "Note that  the `QubitOperator` \n",
    "constructor simplifies `hamil` automatically.\n",
    "Hamiltonians are Hermitian, so after `QubitOperator`\n",
    "finishes simplifying `hamil`, the coefficient of every \"term\"\n",
    "of `hamil` must be real."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hamil=\n",
      " 0.7 [X1 Y2] +\n",
      "0.4 [Y1 Y3]\n"
     ]
    }
   ],
   "source": [
    "hamil = QubitOperator('X1 Y3 X1 Y1', .4) + QubitOperator('Y2 X1', .7)\n",
    "print('hamil=\\n', hamil)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what is the purpose of this Hamiltonian `hamil`?\n",
    "The cost function to be minimized is defined as the mean value of `hamil`.\n",
    "More precisely, if $\\ket{\\psi}$ is the ouput of the circuit\n",
    "specified above, then the cost function equals $\\bra{\\psi} | H | \\ket{\\psi}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial values for parameters `#1` and `#2` to be minimized\n",
    "must also be passed in inside an input dictionary to \n",
    "the constructor of class `MeanHamilMinimizer`. Variable `all_var_nums` contains the keys \n",
    "of that dictionary.\n",
    "\n",
    "> Note: Results are very heavily dependent on initial x_val, probably because,\n",
    "due to the periodic nature of qubit rotations, there are many local minima"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_var_num_to_rads = {1: 1., 2: 3.}\n",
    "all_var_nums = init_var_num_to_rads.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, we define a function `case()` that creates an\n",
    "object of `MeanHamilMinimizer`, asks that object to minimize the cost function,\n",
    "and then `case()` returns the final result of that minimization process.\n",
    "As advertised at the beginning of this notebook, we \n",
    "use various sub-methods of `scipy.optimize.minimize`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_samples = 0\n",
    "print_hiatus = 25\n",
    "verbose = False\n",
    "np.random.seed(1234)\n",
    "\n",
    "emp_mhamil = MeanHamil_native(file_prefix, num_qbits, hamil,\n",
    "            all_var_nums, fun_name_to_fun, simulator_name='SEO_simulator', num_samples=num_samples)\n",
    "targ_mhamil = MeanHamil_native(file_prefix, num_qbits, hamil,\n",
    "            all_var_nums, fun_name_to_fun, simulator_name='SEO_simulator') # zero samples\n",
    "def case(**kwargs):\n",
    "    return MeanHamilMinimizer(emp_mhamil, targ_mhamil,\n",
    "                 all_var_nums, init_var_num_to_rads,\n",
    "                 print_hiatus, verbose).find_min(minlib='scipy', **kwargs) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We focus on two cases, num_samples=0 and num_samples > 0. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When `MeanHamil` is told that \n",
    "num_samples=0, it does no sampling. It just calculates\n",
    "the mean value of `hamil` \"exactly', using the exact\n",
    "final state vector calculated by Qubiter's `SEO_simulator` class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_val~ (#1, #2)\n",
      "iter=0, cost=-0.129256, x_val=1.000000, 3.000000\n",
      "iter=25, cost=-0.128526, x_val=2.275041, 3.209272\n",
      "iter=50, cost=-0.247852, x_val=1.105711, 3.540960\n",
      "iter=75, cost=-0.241892, x_val=0.978839, 3.768145\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   direc: array([[ 0.        ,  1.        ],\n",
       "       [-0.11503823,  0.22534013]])\n",
       "     fun: -0.2479366921073845\n",
       " message: 'Optimization terminated successfully.'\n",
       "    nfev: 98\n",
       "     nit: 4\n",
       "  status: 0\n",
       " success: True\n",
       "       x: array([1.09458962, 3.55198269])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min_method = 'Powell'\n",
    "emp_mhamil.num_samples = 0\n",
    "case(method=min_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When num_samples > 0, Qubiter calculates the final state vector\n",
    "of the circuit, then it samples that num_samples times,\n",
    "obtains an empirical probability distribution from that,\n",
    "and then it calculates the mean value of `hamil` \n",
    "using that empirical distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_val~ (#1, #2)\n",
      "iter=0, cost=-0.099200, x_val=1.000000, 3.000000\n",
      "iter=25, cost=-0.196800, x_val=1.248064, 3.472136\n",
      "iter=50, cost=-0.219200, x_val=1.223266, 3.798374\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   direc: array([[1., 0.],\n",
       "       [0., 1.]])\n",
       "     fun: -0.22179999999999994\n",
       " message: 'Optimization terminated successfully.'\n",
       "    nfev: 61\n",
       "     nit: 2\n",
       "  status: 0\n",
       " success: True\n",
       "       x: array([1.22326557, 3.77708762])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min_method = 'Powell'\n",
    "emp_mhamil.num_samples = 1000\n",
    "case(method=min_method)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_val~ (#1, #2)\n",
      "iter=0, cost=-0.139800, x_val=1.000000, 3.000000\n",
      "iter=25, cost=-0.128000, x_val=0.999909, 2.999462\n",
      "iter=50, cost=-0.158200, x_val=1.009999, 2.999992\n",
      "Warning: Desired error not necessarily achieved due to precision loss.\n",
      "         Current function value: -0.163000\n",
      "         Iterations: 1\n",
      "         Function evaluations: 75\n",
      "         Gradient evaluations: 16\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "      fun: -0.16300000000000003\n",
       " hess_inv: array([[1.06504473, 0.10933946],\n",
       "       [0.10933946, 0.01122658]])\n",
       "      jac: array([ 1.36, -0.4 ])\n",
       "  message: 'Desired error not necessarily achieved due to precision loss.'\n",
       "     nfev: 75\n",
       "      nit: 1\n",
       "     njev: 16\n",
       "   status: 2\n",
       "  success: False\n",
       "        x: array([0.99999857, 2.99999152])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# very sensitive to eps, seems to be hitting discontinuity\n",
    "min_method='CG'\n",
    "num_sample= 0\n",
    "case(options={'disp': True, 'eps':1e-2})"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "12px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
